Sex-specific expression of pheromones and other signals in gravid starfish

Background Many echinoderms form seasonal aggregations prior to spawning. In some fecund species, a spawning event can lead to population outbreaks with detrimental ecosystem impacts. For instance, outbreaks of crown-of-thorns starfish (COTS), a corallivore, can destroy coral reefs. Here, we examine the gene expression in gravid male and female COTS prior to spawning in the wild, to identify genome-encoded factors that may regulate aggregation and spawning. This study is informed by a previously identified exoproteome that attracts conspecifics. To capture the natural gene expression profiles, we isolated RNAs from gravid female and male COTS immediately after they were removed from the Great Barrier Reef. Results Sexually dimorphic gene expression is present in all seven somatic tissues and organs that we surveyed and in the gonads. Approximately 40% of the exoproteome transcripts are differentially expressed between sexes. Males uniquely upregulate an additional 68 secreted factors in their testes. A suite of neuropeptides in sensory organs, coelomocytes and gonads is differentially expressed between sexes, including the relaxin-like gonad-stimulating peptide and gonadotropin-releasing hormones. Female sensory tentacles—chemosensory organs at the distal tips of the starfish arms—uniquely upregulate diverse receptors and signalling molecules, including chemosensory G-protein-coupled receptors and several neuropeptides, including kisspeptin, SALMFamide and orexin. Conclusions Analysis of 103 tissue/organ transcriptomes from 13 wild COTS has revealed genes that are consistently differentially expressed between gravid females and males and that all tissues surveyed are sexually dimorphic at the molecular level. This finding is consistent with female and male COTS using sex-specific pheromones to regulate reproductive aggregations and synchronised spawning events. These pheromones appear to be received primarily by the sensory tentacles, which express a range of receptors and signalling molecules in a sex-specific manner. Furthermore, coelomocytes and gonads differentially express signalling and regulatory factors that control gametogenesis and spawning in other echinoderms. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-022-01491-0.

Pacific COTS (A. planci cf. A. solaris) aggregating in captivity release a complex array of hundreds of proteins and other factors into the water column that attract other individuals, suggesting that pheromones contribute to the formation of spawning aggregations in nature [28]. The reception and processing of these water-borne signals are likely to occur via receptors on external-facing tissues and organs, including the sensory tentacles located at the distal end of each arm [100][101][102]. Here, we use a large-scale transcriptomics approach to identify factors contributing to reproductive aggregations and spawning. We have analysed multiple somatic tissues and organs, and gonads of gravid, wild COTS just prior to a mass spawning event. This approach is in contrast to most previous studies, which have focussed solely on gonads in animals translocated from their natural habitat [103][104][105][106][107][108]. We find that males and females differentially express putative attractants [28] in their gonads and chemosensory GPCRs and other receptors in their sensory tentacles. Females also upregulate a suite of neural signalling molecules in the sensory tentacles. Other neuroendocrine signalling factors involved in regulating reproduction and spawning in other echinoderms [34, 37-41, 46, 47, 49, 50, 103, 104, 108-113] are differentially expressed, including the unexpected upregulation of RGP in coelomocytes [57][58][59][60][61][62][63][64][65][66].

Consistent and distinct gene expression profiles in tissues and organs from wild, gravid COTS
We generated 1.1 billion high-quality reads from 104 CEL-Seq2 transcriptomes from RNA isolated from eight tissues and organs removed from seven females and six males collected from the wild [114]. Tissues and organs were placed in RNALater within 2 h of the individual being dislodged from the reef, minimising transcriptional changes associated with translocation [115]. All COTS had large ripe gonads, with three individuals spawning upon collection (M. Jönsson and M. Morin, personal observation), suggesting that a natural spawning event was imminent [92]. One male radial nerve transcriptome was not further analysed because of a low mapping rate (< 0.5 million mapped reads). The remaining 103 transcriptomes had a 69.3% mapping rate to protein-coding sequences (CDS) in the GBR v1.1 genome (Additional file 1: Table S1.1) [28,115]. A total of 18,032 of the 24,071 CDS were expressed in at least one tissue or organ, with each tissue/organ expressing on average 15,306 CDS (64% of the genome). The majority of all CDS (12,658 CDS; 71% of all expressed CDS) were detected in all tissues/organs.
Principal component analyses (PCAs) and hierarchical clustering of transcriptomes reveal strong clustering associated with tissue/organ type, regardless of sex, with the first principal component (PC1) accounting for 34% of the variation among samples (Fig. 1A). Gonad transcriptomes are distinct from all somatic tissues/organs, and the internally located coelomocytes are transcriptionally distinct from the two sets of external-facing tissues. The two clusters of external-facing tissues separate into oral tissues/organs, which face the substrate (radial nerve, sensory tentacles and tube feet) and aboral tissues/organs, which face the water column (papulae, skin and spines). Separate PCAs of oral and aboral transcriptomes distinguish each tissue/ organ and indicate sensory tentacle, and tube feet transcriptomes are the most similar to each other (Fig. 1B). A global pairwise comparison of all transcriptomes is consistent with the PCAs and distinguishes male and female gonads from each other and the somatic tissues/ organs (Fig. 1C).
The three oral tissues/organs uniquely upregulate 431 CDS (DESeq2 analyses; p-adj < 0.05) (Additional file 1: Table S1.2). These CDS are enriched in 83 GO terms consistent with the neural and chemosensory roles of these organs, such as ion channel and transmembrane signalling receptor activity (Additional file 1: Table S1.3). The three aboral tissues/organs have 382 significantly upregulated CDS (p-adj < 0.05) (Additional file 1: Table S1.4), which are enriched in eight GO terms consistent with these tissues/organs producing extracellular gene products (Additional file 1: Table S1.5). Coelomocytes have 1025 significantly upregulated CDS (Additional file 1: Table S1.6) enriched in 97 GO terms that reflect their diverse roles in metabolism, immunity and cell movement, including proteolysis and cytoskeletal organisation (Additional file 1: Table S1.7). The gonads have 2181 uniquely upregulated CDS (Additional file 1: Table S1.8) enriched in 423 GO terms, many of which are related to cellular metabolism and reproductive processes, including meiosis (Additional file 1: Table S1.9).

COTS tissues and organs have sexually dimorphic gene expression
A comparison of gene expression in gravid male and female COTS using DESeq2 reveals 42% of the expressed CDS (7607) across all tissues/organs are significantly differentially expressed between males and females ( Table 1; Additional file 1: Table S1.10). A majority of these (7448) are differentially expressed between the testes and the ovaries, with 4111 CDS upregulated in the testes. CDS upregulated in the testes and ovaries are significantly enriched in GO terms related to meiosis, and spermatogenesis and oogenesis, respectively. Doublesex and mab-3-related transcription factor 2 (DMRT2), and oral (bottom) tissues/organs. C Pairwise correlation of eight tissue/organ transcriptomes; female and male gonads and coelomocyte transcriptomes are distinct from all other tissues. Pink, tube feet; purple, sensory tentacles; green, radial nerve; turquoise, skin; brown, papulae; blue, spines; red, coelomocytes; yellow, gonads; circles, females; triangles, males. Depicted values illustrate the scaled (z-score) expression levels based on collapsed VST-normalised read counts kisspeptin and bindin are upregulated in the testes, while egg bindin receptors, egg coat matrix proteins and cyclin B are upregulated in the ovaries (Additional file 1: Table S1.10-S1.12).
All somatic tissues/organs surveyed in this study have sexually dimorphic gene expression profiles, with 283 CDS being significantly differentially expressed between sexes, including 14 transcription factors, 16 neuropeptides and five GPCRs (Tables 1 and 2; Additional file 1: Table S1.10). The sensory tentacles have the largest number of differentially expressed CDS (146), including the majority of the differentially expressed transcription factors, neuropeptides and GPCRs (Table 2).

Sex-biassed expression of putative attractants
An aquarium-based study of Pacific COTS previously identified 94 proteins in an exoproteome that is secreted into seawater from non-gravid aggregating Pacific COTS containing conspecific chemoattractants [28]. Here, we find that 84 (89.4%) of the exoproteome CDS are expressed in gravid wild COTS ( Fig. 2A), with 24 and 14 being significantly upregulated in females and males, respectively (Additional file 2: Table S2.1). All but one of these are differentially expressed between the ovaries and testes; beta-d-xylosidase 2 is upregulated in female coelomocytes (Fig. 2B). Among exoproteome CDS differentially expressed in the gonads, five of six differentially expressed ependymins are upregulated in the ovaries. Also, among the known exoproteins, the testes upregulate potential extracellular regulators, otoancorin and transforming growth factor-betainduced protein ig-h3.
In addition to the exoproteome, we identified a further 194 genes encoding secreted factors that are significantly upregulated in male tissues/organs (Additional file 2: Table S2.2), of which 68 are most highly expressed in the testes ( Fig. 2C; Additional file 2: Table S2.3). These include 33 conserved signalling and structural proteins, such as alpha-2-macroglobulin, bindin and beta-microseminoprotein, as well as 15 conserved uncharacterised proteins and two novel proteins.
Putative receptors of conserved neurohormones, including receptors for kisspeptin, gonadotropinreleasing hormone (GnRH), corazonin (CRZ), orexin and RGP are expressed in a wide range of tissues, including the coelomocytes (Fig. 4C). COTS express Depicted values illustrate scaled (z-score) expression levels based on TPM normalised reads from seven females and six males (except radial nerve, which is from five males). Pie charts depict the proportion of highly expressed genes, per tissue, that fall into broad secreted protein classes (blue, hydrolytic enzymes; light green, other enzymes; orange, enzyme inhibitors; green, structural/signalling proteins; red, conserved uncharacterized proteins; yellow, uncharacterised proteins). B Individual COTS expression profiles of the 37 putative attractants that are significantly differentially expressed between the testes and the ovaries. CDS names are listed to the right. Left heatmap, quartile analysis of the mean transcript abundance; Q4, most highly expressed CDS; Q1, lowest expressed CDS; and 0, expression was not detected. Right heatmap, TPM normalised expression. Star (*) indicates the COTS-specific ependymin-related proteins. C Heatmap of 194 putative secreted proteins upregulated in male tissues/organs, 68 of which are highly expressed in the testes. Highlighted are conserved CDS of interest (Additional file 2: Table S2.3) Fig. 3 Neuropeptidome expression in the sensory tentacles and other tissues and organs. A Expression of 39 neuropeptides [122]. Male and female expression profiles are grouped by tissues/organs. Left heatmap, quartile analysis of the mean transcript abundance; Q4, most highly expressed CDS; Q1, lowest expressed CDS; and 0, no expression was detected. The right heatmap illustrates the scaled (z-score) expression levels based on TPM normalised read counts. Star (*) indicates the significant differential expression between male and female sensory tentacles. B Expression level of RGP across the different tissues/organs. Males and females are grouped by tissue/organ. Values depict square root-transformed TPM read counts of RGP. The expression of RGP is significantly higher in the coelomocytes of both sexes compared to all other tissues (p-adj < 0.001). C Volcanoplot of 15,777 expressed CDS in male and female sensory tentacles. Grey dots depict non-significantly differentially expressed CDS, and black dots depict significantly differentially expressed CDS (146) (p-adj < 0.05). Highlighted CDS names: green, neuropeptides; red, transcription factors; purple, GPCRs. CDS with negative and positive Log 2 fold change are upregulated in female and male sensory tentacles, respectively 15 kisspeptin receptors (ApKPRs), of which the majority are highly expressed in the radial nerve. Notably, female sensory tentacles are upregulating the kisspeptin precursor type and a putative cognate receptor (ApKPR11.2) compared to males ( Fig. 4C; Additional file 6: Fig. S2). Putative receptors for CRZ and GnRH are lowly expressed, however appear to concentrate in the spines and radial nerve, respectively ( Fig. 4C; Additional file 7: Fig. S3). Putative receptors for orexin are most highly expressed in the radial nerve, papulae and spines ( Fig. 4C; Additional file 7: Fig. S3). Three putative RGP receptors are expressed in gravid COTS and appear to be echinoderm-specific (Fig. 4C, D). The putative receptors for RGP are lowly expressed in several tissues/organs; however, gbr.152.19.t1 is significantly upregulated in the testes compared to the ovaries (p-adj < 0.05) ( Fig. 4C; Additional file 1: Table S1.10; Additional file 7: Fig. S3).

Procuring wild transcriptomes
In this study, we analysed the transcriptional states of tissues and organs using replicated and defined biopsies from multiple gravid male and female COTS sampled directly on the Great Barrier Reef [114]. To our knowledge, this provides the first comprehensive sex-specific gene expression profiles in somatic tissues/organs and gonads of any echinoderm in the wild. The biologically replicated sampling enables the identification of genes that have consistent differences in the expression between tissues/organs and sexes. The high similarity of tissue transcriptomes between individuals that were collected over 7 days from two adjacent reefs suggests that the sampling regime has accurately captured native gene activity in wild COTS. Importantly, the sampling of COTS directly after their removal from the reef avoids transcriptional changes that may occur during translocation and being held in captivity. We have recently demonstrated that COTS translocated and kept in captivity exhibit large-scale and sustained transcriptional changes, with the expression of over 20% of the CDS in the genome changing significantly from the native expression profiles [115].

Sexual dimorphism in crown-of-thorns starfish
Our analyses of the gene expression in both external and internal somatic tissues/organs, and in the testes and ovaries, provide new insights into the physiological and sensory states of this synchronous broadcast spawner [2]. For instance, we find that males predominantly upregulate secreted conspecific signalling factors in their testes, while females upregulate hormones and neuropeptides in their sensory tentacles. Most of the 94 exoproteins comprising a secretome previously shown to attract COTS [28] are upregulated in specific tissues, including the gonads, where 37 of these putative attractants are significantly differentially expressed between males and females. Moreover, RGP, a known gonadotropin and regulator of starfish reproduction [57][58][59][60][61][62][63][64][65][66], is significantly upregulated in the coelomocytes, suggesting a vital role for these cells in reproduction.
Although most echinoderm species are gonochoristic, there tend to be few visibly discernible external differences between sexes [123][124][125][126]. This is also the case for COTS [127,128]. Nonetheless, males tend to spawn before females [80,85,86], indicating physiological and behavioural differences between the sexes. Our analysis of seven somatic tissues and organs, and gonads, uncovered sex-specific gene expression profiles in all tissues, with the gonads contributing 96.3% of the significant transcriptional differences between males and females. Within the gonads, 46.5% of all expressed CDS are significantly differentially expressed between sexes. The remaining 3.7% of sexually dimorphic gene expression is in somatic tissues and organs, consistent with gravid males and females having distinct physiologies that allow them to differentially respond to environmental and conspecific signals. In particular, the largest sex-specific differences in COTS appear in the most distal organs-sensory tentacles and spines-suggesting that males and females are most different in tissues/ organs that are likely to first receive external molecular signals. There also appear to be sex-specific differences in internal physiology and signalling, indicated by sex-specific differences in coelomocyte gene expression. The differential spawning between sexes in other echinoderms [81,83,84,87,129] suggests that sexual dimorphisms, at least at the molecular level, are more widespread in this phylum than previously appreciated. Differentially expressed CDS in COTS gonads include sex-specific proteins that are conserved in echinoderms. For instance, bindin-an acrosomal protein essential for sea urchin fertilisation-receptors for egg jelly, and guanylate cyclase are all upregulated in the testes [36-39, 109, 110, 113], and egg bindin receptor, egg coat matrix proteins and a putative maturation promoting factor are upregulated in the ovaries [34,111,112,[130][131][132]. In addition, DMRT-like transcription factor, kisspeptin and receptors for melatonin and thyrotropin-releasing hormone are upregulated in the testes. All of these appear to play a conserved role in sex determination and differentiation in deuterostomes and echinoderms [133][134][135][136][137][138].
In aquaria, aggregations of COTS produce complex chemical plumes comprising hundreds of proteins that attract naive solitary adults [28]. Of the 94 characterised proteins in this plume, 84 are expressed in wild COTS, with 38 of these significantly differentially expressed between sexes. This suggests that the attractant proteins identified in the ex situ plumes may be playing a similar role in the wild. Among these are six ependymin-related proteins that are part of a COTS-specific expansion of this rapidly evolving gene family, which belongs to a larger metazoan ependymin family and is related to the developmental gene LVN1.2 in the sea urchin Lytechinus variegatus [28,144,145]. These are strong candidates as conspecific signalling factors. Five ependymin-related proteins are upregulated in the ovaries and one in the testes. In addition to these characterised exoproteins, a large suite of secreted CDS is significantly upregulated in the testes, including a number of conserved factors such as alpha-2-macroglobulin and beta-microseminoprotein. Alpha-2-macroglobulin is a proteinase inhibitor that can bind many hormones [146][147][148], is present in seminal plasma [149] and appears to be involved in conspecific communication in the barnacle Balanas amphitrite [150]. Beta-microseminoprotein also appears to affect mating behaviours in the squid Loligo pealeii [151].
The presence of previously identified proteins in an exoproteome that attracts other COTS [28] and the identification of a raft of new secreted proteins, many of which are highly expressed in the testes, supports the premise that male and female COTS communicate their sexual identity and readiness to spawn through waterborne pheromones. Similar communication systems have been proposed for holothurians [87,152,153], ophiuroids [129], other asteroids [82] and other marine invertebrates [6,151,154,155].

Receptors and neuropeptides may regulate the response to external cues
Pheromones and other signals can affect conspecifics through the activation of receptors and downstream signalling pathways [156][157][158][159]. In wild COTS, we find that the sensory tentacles express the most GPCRs, including class A rhodopsin-like and putative olfactory receptors, consistent with their role as a chemosensory organ [28,[100][101][102]. The two class A rhodopsin-like GPCRs differentially expressed between male and female sensory tentacles are ApKPR11.2 and SSTR2. In addition, the sensory tentacles differentially express a large number of neuropeptides including orexin, which in mammals appears to regulate the sensitivity of expressed olfactory receptors to specific chemical signals [160][161][162][163][164]. Moreover, several of the highly expressed neuropeptides in COTS sensory tentacles and radial nerves have conserved roles in regulating reproduction in echinoderms and vertebrates. For instance, kisspeptin appears to regulate reproductive and metabolic pathways in holothurians [120] and the secretion of gonadotropin-releasing hormone (GnRH) in mammals [116,117]. Mammalian kisspeptin appears to play a role in female-guided attraction towards males [118]. An expanded set of kisspeptin receptors (KPRs) were recently discovered in Asterias rubens and appear to interact with kisspeptin neuropeptides and SALMFamide peptides [165]. A functioning kisspeptin system has also been discovered in the sea cucumber Apostichopus japonicus, which appears to be seasonally expressed, indicating a potential role in reproduction [120]. There are at least 15 putative KPRs expressed in wild COTS. ApKPR1, 3, 8 and 9 are orthologues to KPRs in A. rubens that interact with kisspeptin. ApKPR6 and 7 are orthologues of A. rubens KPRs that interact with SALMFamide, which acts as a muscle relaxant in other starfish [166,167]. The A. rubens orthologue of ApKPR11.2, which is upregulated in the COTS female sensory tentacles, has a low affinity for kisspeptin and SALMFamide [165]. The differential expression of KPRs, kisspeptin and SALMFamide precursors between male and female sensory tentacles suggests that these neuropeptide signalling systems are playing important roles in conspecific communication during reproduction.
COTS coelomocytes appear to be involved in transducing external signals to other tissues and organs, including the gonads. Specifically, we find that the coelomocytes upregulate RGP, which has a conserved role in regulating starfish reproduction [57][58][59][60][61][62][63][64][65][66]. To our knowledge, RGP has not been detected in other echinoderm coelomocytes. The coelomic fluid is a vital system for nutrient transport and immunity in echinoderms [168][169][170] and may also contain neuropeptides and other signals that regulate reproduction and other physiological states [171][172][173][174]. COTS coelomocytes express a large suite of class A rhodopsin-like GPCRs, including a putative corazonin (CRZ) receptor [175,176], KPRs and thyrotropin-releasing hormone receptors. Thus, it appears coelomocytes are competent to receive internally secreted neuropeptides and hormones from the sensory tentacles and radial nerves and are playing a central role in regulating reproduction and spawning.

Proposed regulation of COTS spawning based on gene expression in the wild
The tissue-, organ-and sex-specific gene expression profiles presented in this study were procured from wild, gravid COTS prior to a spawning event. A large number of putative attractant pheromones, spawning-inducing factors and neuropeptides and their cognate receptors are expressed in reproductive and sensory organs, suggesting their involvement in regulating COTS spawning. The testes appear to be the largest source of pheromones that could be released prior to, or at the time of, spawning (Fig. 5). Other putative pheromones appear to originate from the ovary and somatic tissues/organs of both sexes. The sensory tentacles appear to be the primary organ to detect exogenous signals, including pheromones, and to integrate these signals and promulgate information inwards. We posit that a combination of pheromones, receptors and internal signals, including conserved neurohormones and neuropeptides, allows for the detection of factors that guide the formation of aggregations and induce spawning (Fig. 5). Females activate a unique suite of receptors and signals consistent with each sex having a unique response to external cues.
We propose that exogenous signals received by the female sensory tentacles induce the release of kisspeptin, which is differentially upregulated in this organ. This may regulate the secretion of GnRH and/or CRZ from the radial nerve into the coelomic fluid, where they interact with receptors on the coelomocytes and, in turn, induce the release of RGP stores in these cells (Fig. 5). RGP may interact with RGP receptors in the gonads to regulate the last stages of gamete maturation and spawning, as in other asteroids [57][58][59][60][61][62][63][64][65][66][67][68][69][70][71]. To confirm this complete molecular cascade of events will require experimentation ex situ. Our transcriptome dataset on wild gravid COTS provides a framework to guide these future analyses.

Conclusions
Analysis of tissue and organ transcriptomes that reflect the physiological state of gravid COTS around the time of spawning in the field suggests that there is extensive tissue-level sexual dimorphism beyond the gonads and that these differences play a role in COTS aggregation and spawning. Of particular note is the expression of an overabundance of putative pheromones in the testes, including many known to be secreted by aggregating COTS [28], and the expression of a large repertoire of receptors and signalling factors in the female sensory tentacles. These findings suggest mechanisms by which females could detect and respond to male spawning, through a signalling cascade that includes conserved reproductive neuropeptides and that ultimately induces oocyte maturation and female spawning.

Sampling of wild-caught COTS tissues
Seven gravid female and six gravid male COTS were collected from Davies and Lynch's Reefs (18° 50′ S, 147° 39′ E, and 18° 76′ S, 147° 63′ E, respectively) on the Great Barrier Reef between 4 and 9 December 2019 and immediately transferred to aquaria containing ambient seawater onboard a vessel moored at the reef. All tissue and organ samples were dissected and placed in RNALater within 2 h of the COTS being removed from the reef. These were left in RNALater at 4 °C for 24 h and then stored at −20 °C.
Eight different tissues/organs-coelomocytes, gonads, sensory tentacles, tube feet, radial nerves, skin, papulae and spines-were sampled. Coelomocytes were separated from 1.5 ml of fluid extracted from the coelom using a syringe and 21-gauge needle by immediately centrifuging at 2000 × g at ambient temperature. The coelomic fluid was decanted, and coelomocytes were resuspended in RNALater [177]. The gonads were extracted via a small incision at the base of a random arm and gently removed via forceps. The sex of the animals was recorded, and the gonads were placed in RNALater. From each individual, (i) two sensory tentacles were dissected from the tip of a random arm, (ii) one tube foot was dissected from the middle of a random arm, (iii) a 30-mm length of the radial nerve cord was dissected from the middle of a random arm, (iv) the skin was dissected from the proximal base of a random arm, (v) approximately 50 papulae were dissected from the dissected skin samples while in RNALater (the papulae-less skin was used to make the skin libraries) and (vi) the spine samples consisted of the epithelium layer scraped off of one spine, collected from the base of a random arm.

CEL-Seq2 library construction and sequencing
RNA from each individual tissue/organ was extracted using TRI Reagent (Sigma) following the manufacturer's protocol. The quantity and quality of the RNA were assessed using a Qubit ® fluorometer (Invitrogen-Life Technologies) and the Agilent 2100 Bioanalyzer (Agilent Technologies). CEL-Seq2 library construction and sequencing were performed as described in Hashimshony et al. [178] with each individual sample being barcoded. The eight pooled tissue/organ libraries were sequenced on the Illumina HiSeq X ten platform.

Differential gene expression analyses
To remove potential technical read errors often associated with lowly expressed genes [183], we tested different expression threshold levels of an average of ≥ 0.25 and ≥ 1 read per gene per library for all replicates of a given tissue and opted for the expression threshold of ≥ 0.25 mean reads per tissue/organ (Additional file 8: Fig. S4). DESeq2 [184] was performed to identify differentially expressed CDS with a p-adjusted (p-adj) value of < 0.05. Gene Ontology (GO) enrichment analyses of differentially expressed CDS were performed using Fisher's exact test function available on Omicsbox (v1.4.11), using the GO annotation provided by Blast2GO. GO terms were considered enriched with a false discovery rate (FDR) value of < 0.05. Principal component analyses (PCAs) were used to visualise the differences in gene expression between tissues/organs and sexes. PCAs were performed in RStudio on variance-stabilise transformed (VST) counts obtained with DESeq2. Heatmaps were created to visualise the expression patterns of genes of interest, using transcripts per million (TPM) normalised raw reads, in the R package "pheatmap" [185,186]. TPM reads were also used to calculate the quartile values for each individual. Expressed genes were classified into four quartiles ranging from lowest (Q1; bottom 25%) to highest (Q4; top 25%) level of transcript abundance. Upset plots were used as an alternative to Venn diagrams, to illustrate the eight-way expression of GPCRs between the eight different tissues. These upset plots were constructed using the R package UpsetR [187]. All analyses and visualisations were performed in RStudio version 4.0.2 [188].

Identification of secreted factors
The COTS v1.1 gene models (GBR v1.1) [115] were screened for genes encoding signal peptides using Sig-nalP 5.0 [189] and then for transmembrane domains using TMHMM Server v. 2.0 [190,191]. Proteins were considered secreted if they possessed a predicted signal peptide but no transmembrane domain. CDS were blasted on NCBI using blastp, and sequences with no significant matches (E-value < 0.05) were considered potentially novel species-specific proteins.

Additional file 5:
Additional file 8: Fig. S4. Comparison of different expression threshold values. A PCA of all individual transcriptomes separated into respective tissue cluster (colour) and sex (shape). The threshold value of expression was set to a mean number of ≥ 0.25 reads per tissue. B Same PCA as (A) but here the threshold value was set to a mean number of ≥ 1 reads per tissue. This comparison reveals a small difference between the two threshold values. (A) resulted in 18,032 CDS genes counted as expressed and (B) resulted in 16,127 CDS genes counted as expressed. Due to the replicated sampling, we opted to include as many CDS as possible in the differential expression analysis.
Additional file 9: Fasta file of sequences used to infer RGP receptor phylogeny.